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We numerically investigate a particular type of spin solitons inside a trapped atomic spin-1 Bose- 
Einstein condensate (BEC) with ferromagnetic interactions. Within the mean field theory approxi- 
mation, our study of the solitonic dynamics shows that the solitonic wave function, its center of mass 
motion, and the local spin evolutions are stable and are intimately related to the domain structures 
studied recently in spin-1 *^Rb condensates. We discuss a rotating reference frame wherein the 
dynamics of the solitonic local spatial spin distribution become time independent. 
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I. INTRODUCTION 

Solitonic structure is an interesting and common sub- 
ject in nonlinear science. It has been extensively studied 
in nonlinear optics 0, fluids 0] and plasmas mag- 
netic films 0, nanotubes and recently in quantum 
superfluid of atomic Bose-Einstein condensates |^. In 
nontechnical terms, the shape of a soliton remains unal- 
tered during propagation due to a balance of the spread- 
ing from its dispersion with nonlinear interactions. This 
feature has been applied with great success in optical 
soliton based advanced communications. 

As realized in recent years an atomic Bose-Einstein 
condensate (BEC) provides a new and powerful tool for 
the study of nonlinear phenomena. At weak interac- 
tions, a condensate is successfully modelled by the Gross- 
Pitaevskii (GP) equation, or the nonlinear Schrodinger 
equation (NLSE) S B E El II3- Feshbach reso- 
nance in recent years has provided a practical mechanism 
for adjusting the atomic s-wave scattering length through 
which both the sign and the strength of the atomic in- 
teractions become adjustable Bright matter wave 
soliton trains were realized in a ^Li condensate due to its 
attractive atom-atom interaction, while gray and dark 
solitons were realized for repulsively interacting conden- 
sates through the phase en gine ering technique on the con- 
densate wave function 

[Hill. 

Solitons with more than one parameter are often re- 
ferred to as vector solitons. For instance, vector solitons 
are created in nonlinear optics if different polarizations 
or spatial modes of light in an optical fiber are taken 
as parameters. This could easily be implemented in a 
spinor condensate as the atomic internal degrees of free- 
dom become available. Many have suggested creating 
a Manakov soliton in a two-component BEC with at- 
tractive self-interactions Symbiotic solitons have 
also been investigated in a two-component BEC with re- 
pulsive self-interactions accompanied by attractive cross- 
interactions E3- Domain wall solitons can be created 
in a trapped two-component condensate even when both 
self- and cross-interactions are repulsive More gen- 
erally, self-consistent collective excitations of an atomic 



condensate, or higher eigen-states of the Gross-Pitaevskii 
equation, are also called solitons loosely because of their 
reasonable stabilit y a. nd resemblance to other properties 
of solitonic states [l9l |. 

With a far off resonant optical trap, all hypcrfinc states 
of an atom can be confined simultaneously, thus realizing 
a spinor condensate, as has been reported in many ex- 
periments with 23Na and *'^Rb atoms "l^, IH Hi IH 
|24, ^ ^ m il. The releasing of the spin degree 
of freedom has led to many novel properties lacking in 
a single component or scalar condensate. For instance, 
the spin exchange interaction as in a spin-1 condensate, 
whereby two \F — l,mF — 0) atoms are converted into 
one \F = 1, mp — +1) atom, and one \F — 1, mp — —1) 
atom upon an elastic collision, already well understood 
|29l ISOl l3ll | , also helps to support the existence of spinor 
solitons. In a homogeneous system, exact soliton solu- 
tions for an attractively interacting spin-1 condensate 
have been recently obtained analytically [32|. Also the- 
oretical investigations have been carried out on mag- 
netic (or spinor) solitons in a ferromagnetically interact- 
ing spin-1 condensate [s^lls^l. 

In this paper, we focus on the more realistic case of spin 
solitons inside trapped atomic condensates (of *^Rb). 
First we attempt to numerically find solitonic structures 
through the propagation in time of the GP equation for 
a trapped spin-1 condensate. Then we demonstrate that 
the dynamics of the found spin solitons is stable and is 
related to the multi-domain structure as studied recently 
|35l Is^ . Finally, we discuss a rotating reference frame 
wherein the local spin distribution of the spinor conden- 
sate becomes stationary like a spatially propagating soli- 
ton in a moving reference frame. 

This paper is organized as follows. In Sec. Uwe give 
an introduction and discuss the motivations for our work. 
The theoretical formulation is based on a model system 
of a spin-1 condensate as is presented in Sec. |n| Sec- 
tion ^O] describes the numerically obtained soliton state, 
while Sec. IIVI focuses on the investigation of the soli- 
tonic state's dynamics, center of the mass motion, and 
local spin distribution. In Sec. we discuss a rotating 
reference frame and show that the local spin distribution 
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becomes time-independent within this frame. Finally we 
conclude in Sec. IVll 



II. REVIEW OF AN EFFECTIVE 
QUASI-ONE-DIMENSIONAL DESCRIPTION 

A gas of interacting spin-1 atoms is described by the 
following second quantized Hamiltonian (summation over 
repeated indices is assumed) [sO, |^ 
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-^V2 + yext(r))vl/, 
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+ ^^l^]^,^. + f *t vl,t (^^)^^. (^^)^^ ^^.^^J ^ (1) 

where ^j{r) (^]) is the field operator that annihilates 
(creates) an atom in the jth internal state at location r, 
j — +, 0, — denotes atomic hyperfine state \F = l,mF = 
+1, 0, —1), respectively. M is the mass of each atom, and 
Kixt(r) is an internal state independent trap potential. 
Terms with coefficients cq and C2 of Eq. (QJ describe 
elastic collisions of two spin-1 atoms expressed in terms of 
the scattering lengths aq (02) in the combined symmetric 
channel of total spin (2), cq = 47r?i^(ao-l-2a2)/3M, and 

=x,y,z are spin-1 matrices. 
The equation of motion for the field operator in the 
Heisenberg picture is given by 



ot 



(2) 



Adopting a mean field theory approach by assuming that 
the condensate consists of a large number of atoms, we in- 
troduce the condensate order parameter or wave function 
= (^i) for the ith component. Neglecting quantum 
fluctuations we arrive at the coupled Gross-Pitaevskii 
equations 
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where £ = ~{fC /2M)\7'^ + V^^tir) + Cf^n, n ^ J2i'^i is 
the total condensate density, and — Another 
way to derive the above coupled GP Eqs. Q is to take 
the variations of the energy functional with respect to 
the condensate wave function 
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Inside a trap with tight radial confinement, the con- 
densate assumes a prolate shape. Several effective one- 
dimensional (ID) approaches have already been devel- 
oped m m m IS El 113: with the simplest of 
them assuming a fixed transverse Gaussian profile. Re- 
cent studies, however, have indicated that the effective 
quasi-lD non-polynomial Schrodinger equation (NPSE) 
is the more appropriate choice. In fact for most sys- 
tems of interest, NPSE represents the most powerful and 
efficient tool especially in the weakly interacting limit 
[3^ m, 113, 113 • In several recent experiments, a single 
running wave optical trap is used to confine spin-1 atomic 
condensates, a situation well described by the quasi-one- 
dimensional trapping geometry [2^. I25I |2^ UM- More 
precisely the external cigar-shaped trap is described by 
a harmonic trap potential Vcxt = {M/2){uj'^r'^ + ujIz^) 
with ujx ~ ijJy = o;^, r±_ = + u'^, and ujz ^ 

The quasi-one-dimensional NPSE description assumes a 
factorized wave function 



N(b^{fr,xiz,t))Mz,t), (5) 



into the transversal and longitudinal functions p2 | , which 
in the case of a spin-1 condensate are governed respec- 
tively, by 
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In the above, N is the total number of condensed atoms. 
(j)i is the factorized longitudinal envelop function of the 
quasi-one-dimensional wave function, and it depends only 
on z and t. <j)±_ is the transversal wave function, satisfy- 
ing J — 1, and is assumed to be identical for 
all three spin components, x is a variational functional 
characterizing the width of the transversal wave func- 
tion, ho = -{h^/2M){d'^/dz^) + V{z) + + coNrjp 
with V{z) — (M/2)w|z^, and E± is the transverse mode 
energy functional. 
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vix) — J df±[(f>±['^ is a scaling factor for the nonlinear 
interaction strength. p{z) — and 6*2 is indepen- 

dent of X and is shown to be given by 
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In obtaining the above relations self-consistently, we have 
also assumed a weak time and z dependence of the trans- 
verse wave function, i.e., d(j>±/dt ~ and V^0^ ~ 
V2 6,. 
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For a condensate with a large number of atoms, as 
in most current experiments, the density distribution of 
the transversal direction approaches the Thomas-Fermi 
(TF) limit, /i ^ ?iw_L, thus we take the TF ansatz for the 
transverse wave function, 



0, Ty_ > X- 



The kinetic energy in the transverse direction is therefore 
neglected under the TF limit, leading to the transverse 
mode energy and scaling factor being 

4 



III. 



SOLITONS IN A TRAPPED SPIN-1 BOSE 
CONDENSATE 



Condensate solitons we discuss in this paper are sim- 
ply self-consistent high energy eigen-states of the GP Eq. 

Armed with the quasi-lD NPSE approach and 
the imaginary time propagation method, we devise the 
following operational procedure for an extensive search of 
the solitonic states. Numerically we first find the ground 
state, then we propagate the quasi-lD NPSE Eq. © in 
the imaginary time domain [42j, continuously projecting 
out the ground state at each time step during evolution. 



|0(")) ^ (0(0) 10(0)), 
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hoping for a self-consistent eigen-solution to be found. It 
is well known that different eigen-solutions of a nonlin- 
ear Hamiltonian may not be orthogonal, furthermore the 
family of all eigen-solutions does not generally constitutes 
an orthornormal and complete basis. Yet, our procedure 
is found to converge at least for the first self-consistent 
collective excitation. This solitonic state we find is essen- 
tially independent of the choice for the transverse mode 
as being a fixed Gaussian or a fixed TF ansatz. 

An example of the spin-1 soliton found in a trapped 
^"^Rb condensate is illustrated in Fig. ^ The scatter- 
ing lengths used are oq = 101.8 and 02 = 100.4 as 
with as being the Bohr radius 44]. In the subplot Fig. 
^a) we display the spatial distribution of the local spin 
average that is defined as 



= {F,)^V2Re[M<P*++r-)]/p, 
= {F.) = {\cb+\'-\cb-\')/p, 



(10) 



We 



whose magnitude is given by f{z) — ^ + fy + fz 
find that the component fy, which is displayed perpen- 
dicular to the plane, remains very close to zero for the 
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FIG. 1: (Color online) The spatial distribution of the local 
spin (a), the density (b), and the phase (c) of a soliton state 
in a trapped spin-1 ^'^Rb condensate. The trap parameters 
are uj^ = ujy ^ (27r)240 Hz, and W2 = (27r)24 Hz. The to- 
tal number of atoms in the trap is N = 10®, and the total 
magnetization {Ai = A'^+ — A'^-) is zero. The solid curves in 
panels (b) and (c) denote the density and the phase for the \+} 
component, the dash-dotted curves denote the |0) component, 
the dashed curves denote the | — ) component, and the dotted 
curve in (b) denotes the total density. The solid line in panel 
(d) shows the local spin angle distribution and the dashed line 
is that for a magnetic kink soliton in a homogeneous system. 



soliton, i.e., the local spins stay almost in the x-z plane. 
We believe the residue small but nonzero value of fy is 
due to the numerical precision. Instead of getting rid of 
these annoying numerical errors, we keep them to study 
the dynamical stability of the soliton state. Figure QJa) 
clearly shows that the local spin rotates nearly a com- 
plete revolution, or 27r around the y-axis from one side 
to the other side of the condensate along the quasi-lD ax- 
ial direction. It is easy to check that the winding number 
of this solitonic spin distribution is about one along the 
y-axis, which shows that the soliton state we find might 
be related to a quasi-one-dimensional vortex. 

We illustrate also the density and the phase distribu- 
tion in Fig. n^b) and^c). We see a n phase shift in 
the |0) component near the central region where the den- 
sity of the |0) component diminishes. Also we notice 
that the |-|-) and |— ) components are symmetric with 
respect to the trap center at z = 0. The panel (d) 
shows the angular distribution of the local spin (solid 
line), cos6' = fz/f- It is compared to the soliton angu- 
lar distribution (dashed line) of a magnetic kink soliton 
given by = 2 arctan(z/^) where ^ = 5.40^ and the phase 
of |0) component jumps in its neighborhood. The good 
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FIG. 2: (Color online) Temporal propagation of the spatial 
profile of the trapped soliton. The parameters and notations 
are the same as in Fig. The upper and the middle rows dis- 
play the density distributions of the soliton at different times. 
The lower row shows the time dependence of the fractional 
populations of every spin components. 



agreement indicates that the soliton state we find in a 
trapped spin-1 condensate is closely related to a mag- 
netic kink soliton. 



IV. REAL TIME PROPAGATION OF THE 
SOLITON IN A TRAP 

Starting from the spinor soliton state found above, we 
can propagate in real time the quasi-lD NPSE in the 
trap. This allows us to easily check the dynamic stabil- 
ity of the solitonic structure. In addition, we can check 
several conservative quantities such as the energy, the 
total number of atoms N , and the total magnetization 
AA = A^_(_ — to confirm the accuracy of our numerical 
integration. During the propagation time t € [0, 5000] 
(l/wz), the relative fluctuations of the energy, the num- 
ber, and the total magnetization are found to be less than 
10"^ 1G"^ and IQ-^^, respectively. 

We show in Fig. |21 the time dependence of the frac- 
tional population for each component, clearly periodic 
in time. The period is found to be T ~ 313 (l/cj^), or 
equivalent to about 2.1 s in real time unit. For a homo- 
geneous condensate this period approaches infinity in the 
thermal dynamic limit, a topic that is beyond the scope 
of the present article and also currently unreachable in 
experiments. For a finite sized condensate as discussed 
here, the period is found to decrease. Figure[21also shows 
typical density distribution for each component at vari- 
ous times during the first half period. We see the |-|-) 
and |— ) components tunnel through each other with the 
assistance of the |0) component and return to their orig- 
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FIG. 3; (Color online) The velocity and the position for the 
center of mass of each soliton component propagating in a 
trap. The parameters and notations are the same as in Fig. 
The motion is periodic, and no noticeable distortion occurs 
after more than 15 periods or about 33 s. 
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FIG. 4: (Color online) An eye shaped diagram for the phase 
portrait of the spinor soliton. The parameters and notations 
are the same as in Fig. Note that the velocity of the |0) 
component is very close to zero, and there are many small 
oscillations for the |-|-) and | — ) components. 



inal locations within each period. We note that the |+) 
and |— ) components are immiscible and the tunnelling is 
generally prohibited without the presence of the |0) com- 
ponent. We do not see such an immiscibility here, since 
the |0) component provides a cohesion for the |-|-) and the 
|— ) components. This type of tunnelling process can also 
be explained in terms of a mutual precession among the 
three components as in Ref. j^^l ■ We actually extract the 
precession angles for both the |-|-) and |— ) components 
and both are found equal to tt, corresponding to a simple 
exchange of their positions after the tunnelling, i.e., spin 
up become spin down and vice versa. The interesting 
feature for a trapped system, distinct from the homoge- 
neous case, is that the tunnelling process repeats itself 
cyclically upon reflections from the trap boundaries. 

Figures 01 and 01 display the position and the velocity 
of the center of mass (COM) for each component of the 
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solitonic condensate Complementary to Fig. [3 Fig. 
|31 shows the oscihations for the last period of our simula- 
tion, which lasts from t = to t — 5000 (l/w^) with the 
oscillation period being about 313 (l/wz). The velocities 
of the COM for the |+) and the |— ) components exhibit 
many small oscillation structures while that of the |0) 
component is always very close to zero. We note from 
Fig. O that the position of the COM for the |0) com- 
ponent is always zero, but the position for the |-|-) and 
the |— ) components oscillates back and forth periodically 
and is symmetrically located on the different sides of that 
of the |0) component. Figure 01 illustrates an eye shaped 
diagram for the phase portrait of the |-|-) and the |— ) 
components during the first period. The "eye" opens 
widely with a clear lid, and the |0) component sits at the 
center. We note that the fluctuations of the velocities 
for both the |+) and the |— ) components are small near 
the center of the trap and become enhanced around the 
turning points of trap boundaries. 

Now we comment on the issue of dynamical stability of 
a ferromagnetically interacting spin-1 Bose condensate. 
Earlier we studied in Ref. 35], the dynamics of a uni- 
formly distributed spin-1 condensate with ferromagnetic 
interactions is unstable. While in this study we see lit- 
tle distortion of the wave function after an extended real 
time propagation even in the presence of initial numerical 
noise. This clearly indicates the soliton state is dynam- 
ically stable. If we compare the soliton state with the 
final domain structure developed from the homogeneous 
spin-1 condensate, we see some similarities to their den- 
sity distributions. These similarities point to a possible 
connection of the trapped solitonic structure we discuss 
to the domain structure in the homogeneous case. On a 
more general term, domain structure can be interpreted 
as a certain restricted kind of multiple solitons. Previ- 
ously, Kasamatsu and Tsubota also noticed these simi- 
larities between the domains and the solitons in a system 
of two-component atomic condensates [4^ . 

We note that the oscillation period is about 2 seconds, 
which is of the order of the experimental lifetime of the 
condensate 0. Thus, it remains challenging to observe 
the spin soliton structure we study here. We also observe 
that the quantum diffusion time of the spin in the ground 
state condensate is only a few trap cycles, usually shorter 
than period of the soliton oscillations [43 . Furthermore, 
our investigation focuses on excited state solitonic struc- 
tures where the application of the single spatial mode 
approximation may become questionable [43 ■ A more 
general treatment that potentially addresses these ques- 
tions seems completely out of reach for the time being. 
On the other hand, as we reduced to iV = 10^, still a 
large number of trapped atoms, the period of the soli- 
ton oscillation reduces due to smaller overlap for the |±) 
components to tunnel through. Thus some aspects of 
these solitonic structures could become observable with 
the system parameters being tuned to optimize the in- 
tended signals before the condensate disappears. 

The detailed local spin evolution of the soliton state is 
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FIG. 5: (Color online) Spatial and temporal distributions of 
the spin evolution for the spinor condensate soliton. The mag- 
nitude of the local spin average is always close to unity due 
to the ferromagnetic interaction, but the local spin direction 
evolves with time due to the spin exchange interaction, (a) 
shows the time dependence of the spatial distribution for the 
azimuthal spin angle i/J and (b) for the polar spin angle 9, 
(c) shows the typical spatial distribution of the spin average 
(t=80), and (d) shows the projection of the local spin distri- 
butions onto the x-y plane. 



further illustrated in Fig. [S| The magnitude of the local 
spin / = ^ fx + fy + fz is nearly the same everywhere 
in space and time. We believe / is exactly the same if a 
more accurate initial state is used. The spatial distribu- 
tion of the pointing angles for the local spin, either the 
azimuthal one ip defined by tanip = fy/fx in|5fa) or the 
polar one 9 defined by cos 6* = /z// inEIb), are strongly 
periodic. Both angles vary rapidly in the central region 
and slowly around the boundaries at every moment (the 
horizontal directions in Fig. Eta) and (b)). Near the 
edge of the condensate (the vertical direction), the angle 
tp essentially remains the same at different times (note 
— TT and TT represent the same spatial direction), while 
the polar angle 6 changes slightly during each period. 
9 changes most rapidly with respect to time at about 
z = ±6a2 where the |-|-) and the |— ) components take 
their highest densities at i = 0. In the central region, 
both angles change moderately with respect to time. 

To further clarify the above analysis for the averaged 
spin directions, we present the 3D spin distribution at 
time i = and t = 80 {1/ujz) in Fig. [Ha) and Fig. 
\^c) and (d) , respectively. As stated previously, the local 
spins have nearly zero projections in the y-axis at time 
t ~ 0, T/2, T, • • •, when they are basically lying in the 
x-z plane and winding around the y-axis with a winding 
number of approximately equal to one [Fig. H^a)]. Simi- 
larly, ~ at time t ~ T/4, 3T/4, 5T/4, • • •, when the 
local spins are lying in the x-y plane but winding around 
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FIG. 7: (Color online) Azimuthal angle of the soliton after the 
transformation. The spatial pattern of the azimuthal angle is 
time independent in a rotating reference frame. 



FIG. 6: (Color online) Density and phase distributions of 
the soliton after the transformation (Eq. IIH . The density 
becomes time independent. 



the z-axis with a similar winding number [Fig. I^c) and 
(d)]. We find the winding number is a topological quan- 
tity that remains conserved during the propagation. 



V. STATIONARY SOLITONS IN A ROTATING 
REFERENCE FRAME 

A propagating sohton with a constant speed is actu- 
ally stationary in a special moving frame. A wonderful 
story about this is that J. S. Russell once followed the 
water wave soliton for a couple of miles. Such impressive 
stability during propagation is also manifested for a soli- 
ton in a spin-1 condensate. We find that the density of 
the soliton becomes time independent with the following 
canonical transformation, 



1 

1 



(11) 



V2 



1 

V2 



^(0+ + (^_)- 



Figure |H| shows the density and phase distributions of 
components |a), |&), and |c), which are all time indepen- 
dent after the transformation. The spatial pattern of the 
phase for each component is also time independent, i.e., 
9a and 9c are even functions along the axial coordinate 
with respect to the center of the condensate, and 9b has 
a TT phase shift at the center of the condensate. We note 
this density distribution is similar to that of a Mermin- 
Ho vortex state of Fig. 1 in Ref. again a potential 
connection of the solitons and the vortex state. 

We can choose a new quantization axis z', different 
from the axial direction of the condensate, such that ipa, 



(j)b, and (f)c correspond to the wave functions of the |-|-), 
|0), and |— ) components in the new coordinate system. 
We then find the spatial distribution oi fz' = I'/'ap — I'/'cP 
is time independent (Fig. EJ, and the fx' and fy' rotates 
with a constant angular velocity around the z' axis as 
shown in Fig. [T] If we go further into the rotating frame 
defined by 



x" ~ x' cosftt + y' sinflt, 
y" ~ — x' sin^t + y' cos ilt, 
z" = z'. 



(12) 



where 51 = 2tt/T, wc find the spatial distribution of the 
spin of the condensate, which is similar to Fig. Q^a), is 
completely time independent. Thus we show again in the 
rotating frame (the double primed one) that the appro- 
priate winding number is time independent and remains 
conserved with change of the reference frame as we con- 
firm previously at the end of the Sec. IIVI 

Further examination shows that the unitary transfor- 
mation matrix of Eq. Hll|l is easily decomposed into spin 
rotations 



U = 




exjp{in) exp(— iFj, • — ) exp(— iF^ • tt). (13) 



The new coordinate {x',y',z') after the rotation cor- 
responds to the original {z,—y,x). Thus fx' = fz, 
fy' = —fy, and fz' = fx- We have checked and con- 
firmed in the original coordinate that fx is indeed time- 
independent, i.e., fz' does not change with time (see Fig. 
EJ. We have also found that the time dependence of 
ta,n~^{—fy/fz) = taTr^[fy' I fx') = cf) \s essentially linear 
as in Fig. [7| Thus we reach the following simplified pic- 
ture as Fig. |S1 in the original coordinate with the given 
initial conditions as shown in Fig. ^a): local spins rotate 
around the x-axis with the same angular frequency Itx jT . 
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FIG. 8: Schematic diagram of the dynamics of local spins in 
the new coordinate system. 



This rotation is constrained so that the x-component of 
the local spin instead of the usual z-component is con- 
served in the discussed soliton state. An amazing aspect 
of this local spin processing dynamics is that the equiv- 
alent gyro-frequency is independent of space. 



VI. CONCLUSION 



tion of the NPSE inside the trapped environment. The 
tunnelling process of the |-|-) and the |— ) components 
with the assistance of the |0) component is explained 
as mutual precessions among the different components. 
We further show that the soliton state is dynamically 
stable and is related to the domain structure observed 
recently in spin-1 ^^Rb condensates j28L l35l|. By trans- 
forming to a rotating reference frame, we demonstrate 
that the time-dependent dynamics of the local spin of 
the soliton state become time-independent. In all of the 
reported numerical simulations we use N = 10^ atoms 
(although unreported simulations were also performed for 
N — 10^ as further checks). Within the mean field the- 
ory we find that the results are basically the same except 
that the period becomes shorter for smaller N due to 
the reduced overlapped region for the |±) components to 
tunnel through. 



In conclusion, we search numerically for a collectively 
excited, or a soliton state in a trapped spin-1 ^''Rb con- 
densate with ferromagnetic interactions by employing the 
quasi-one-dimensional nonpolynomial Schrodinger equa- 
tion (NPSE). We investigate the dynamics of the density, 
the center of mass motion, and the local spin distribution 
of the solitonic state we find based on a real time propaga- 
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